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ABSTRACT 

The addition of metals to any gas can significantly alter its evolution by increasing 
the rate of radiative cooling. In star-forming environments, enhanced cooling can po- 
tentially lead to fragmentation and the formation of low-mass stars, where metal-free 
gas-clouds have been shown not to fragment. Adding metal cooling to numerical sim- 
ulations has traditionally required a choice between speed and accuracy. We introduce 
a method that uses the sophisticated chemical network of the photoionization soft- 
ware, Cloudy, to include radiative cooling from a complete set of metals up to atomic 
number 30 (Zn) that can be used with large-scale three-dimensional hydrodynamic 
simulations. Our method is valid over an extremely large temperature range (10 K ^ 
T ^ 10^ K), up to hydrogen number densities of 10^^ cm~'^. At this density, a sphere of 
1 Mq has a radius of roughly 40 AU. We implement our method in the adaptive mesh 
refinement (AMR) hydrodynamic/N-body code, Enzo. Using cooling rates generated 
with this method, we study the physical conditions that led to the transition from 
Population III to Population II star formation. While C, O, Fe, and Si have been pre- 
viously shown to make the strongest contribution to the cooling in low-metallicity gas, 
we find that up to 40% of the metal cooling comes from fine-structure emission by S, 
when solar abundance patterns are present. At metallicities, Z ^ 10""' Zq, regions of 
density and temperature exist where gas is both thermally unstable and has a cooling 
time less than its dynamical time. We identify these doubly unstable regions as the 
most inducive to fragmentation. At high redshifts, the cosmic microwave background 
inhibits efficient cooling at low temperatures and, thus, reduces the size of the doubly 
unstable regions, making fragmentation more difficult. 

Key words: stars: formation methods: numerical 



1 INTRODUCTION 

The first luminous objects in the universe formed from pri- 
mordial gas, comprised solely of H and He, with only trace 
amounts of D an Li. The relatively simple chemistry of 
metal-free gas, c ombined with tighly constrained cosmolog- 
ical parameters (|Spergel et al.|[2007^ . has allowed the for- 
mation of the first stars to be simulated with extremely 
high precision, from the hierarchical growth of their host 
dark matter halos through to the point where the dense 
proto-stellar cores b e comes optically thick ([Abel et al.ll2002l: 
Bromm etHI |2002|: iBromm fc Loebl 120041: lYoshida et all 



20061 : lO'Shea fc Normanll2007l : iGao et al.ll2007h . With the 
deaths of these stars came the creation of the first heavy 
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elements. Core-collapse and pair-insta bility supe rnovae cre- 
ated metals in copious amount s (Heger fc Wooslev. ,20o3 ) 
and ejected them into the IGM (|Madau et al.ll200ll ). 

The presence of metals alters the dynamics of collaps- 
ing gas-clouds by increasing the number of available atomic 
and molecular transitions, allowing the gas to lose its inter- 
nal energ y more qu ickly tha n in case of no metals (|Omukail 
[2000; Br omm et al.i 2001: B romm fc Loebl l2003bl ). The in- 
troduction of metals adds a new level of complexity to the 
problem of si mulating the forma tion and evolution of cos- 
mic structure. lAbel et all l| 19971 ) identified a minimal set 
of 21 chemical reactions necessary for accurately following 
the non-equilibrium evolution of a g as consisting s o lely o f 
species of H and He, including H2. ICalli fc Fallal (Il998l ) 
showed that 33 total reaction s were required when including 
D and Li species to the gas. lOmukail l|2000l ) performed one 
of the first numerical studies of collapsing gas-clouds to con- 
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sider the contribution of metals. Their chemical network of 
H, He, C, and O included 50 atomic and molecular species 
and 478 reactions. While theirs was not a minimal model, 
the above examples illustrate the great expense associated 
with the expansion of chemical networks to include addi- 
tional elements. Other works have studied the effect of met- 
als on_stax-fomiing gas using similar methodo l ogies to tha t 
of lOmukail hOOdj). e .g.. ISchneider et"all l|2002l . |2003| . l2006h : 
lOmukai et al.l ([20051 ). The complexity of the chemical net- 
works used in these studies limited their treatment of gas 
evolution to one-zone, semi-analytical models. In the earli- 
est work to incorporate metal cooling into three-dimensional 
hydrod ynamic simulat i ons to study metal-enriched star for- 
mation, iBromm et al.l (|200ll ) used a small set of the most 
dominant ato mic transitions of C , N, O, Fe, Si, and S, 
as decribed by iRicotti et al.l l|l997h . Their method also ig- 
nored the cooling from H2, which was justified within their 
study by the assumption of a very large photo-dissociating 
UV background, but is otherwise an extremely important 
coolant in lo w-metallicity environments. For high temper- 
ature gases, [Sutherland Dopital l|l993l ) computed metal 
cooling functions that included 14 heavy elements over a 
range of metallicities, with solar abundance patterns. These 
cooling functions are useful for simulating the IGM and 
other hot, ionized environments, but a minimum tempera- 
ture of 10* K makes them inapplicable to studies of the cold, 
neutral gas associated with star-formation. These cooling 
functions assume coUisional equilibrium of the species and 
as such cannot capture the important role of UV and X-ray 
radiation. 

We introduce a new method for including the cool- 
ing from heavy elements in large-scale hydrodynamic sim- 
ulations that is valid over a wide range of physical con- 
ditions, covers a great number of elemental species, and 
is fast enough to be used in large-scale numerical sim- 
ulations. We have util ized the established photoioniza- 
tion software. Cloudy (|Ferland et al.l Il998h to construct 
large grids of metal cooling data. We have developed a 
method to include both the cooling from heavy elements 
and the non-equilibrium cooling from H2 in hydrodynamic 
simulations. This method has been used successfully in 
the numerical simul a tions of star formation performed by 
ISmith fc SigurdssonI (|2007l ). In ij2] we describe our method 
for creating the metal cooling data, including a new code to 
expedite the process. We, then, present two implementations 
of the cooli ng method in the AMR, hydrodynamic /N-b ody 
code, Enzo (|Brvan fc Normanlll997l : lO'Shea et al.ll2004r ). In 
331 we focus on the application of metals to low-temperature 
environments, identifying the dominant cooling mechanisms, 
and studying the possibility of fragmentation and thermal 
instability in metal-enriched gas. Finally, we end with a dis- 
cussion in 34] of the role played by the heavy elements in the 
formation of structure in the early universe. 



2 NUMERICAL METHOD 

2.1 Calculation of Metal Cooling Rates 

At the current time, it is still too computationally expen- 
sive and memory intensive to follow the non-equilibrium 
chemistry for a large set of heavy elements in a three- 
dimensional hydrodynamic simulation. The exact mass of 



the first massive stars is not known l|Abel et al.l I2OO2I : 
iTan fc McKe3l2004l : lYoshida et al.ll2006h. Also unknown are 
the exact yields of e a rly supernova e ("Hcger & Wooslev"200a; 
Maeder et al] l2005l : llNfomoto et"al, 2006. ; Rockefeller et al] 



20061 ). Similarlv. in many astrophysical systems one might 



want to model computationally the exact metal distribu- 
tions. Consequently, it is not clear a priori what level of 
sophistication of cooling model is needed to adequately cap- 
ture the hydro and thermodynamic evolution of the gas un- 
der consideration. Note that uncertain grain physics also 
increases the potentially important parameter space. In our 
approach, we assume ionization equilibrium, which allows us 
to calculate, in advance, the cooling rate for a parcel of gas 
with a given density and temperature, with incident radia- 
tion of known spectral shape and intensity. Fo r this problem, 
we find the photoionization code. Cloudy (|Ferland et al.l 
1998), especially apt. Cloudy is conventionally used to model 
the transmitted spectrum from a cloud of gas with a given 
chemical composition, being irradiated by a specified source. 
The code must calculate an equilibrium solution by balanc- 
ing the incident heating with the radiative cooling from a 
full complement of atomic and molecular transitions, as well 
as continuum emission from dust. The chemical network of 
Cloudy covers all atomic species from H to Zn, as well as a 
multitude of molecular species. Each elemental abundance 
can be specified individually, giving us the ability to model 
the cooling from a gas with any composition. Since Cloudy 
permits the use of virtually any input spectrum, we are able 
to create cooling data that is suitable for any radiation en- 
vironment. Instead of allowing the code to cycle through 
temperatures until converging on a thermodynamic equi- 
librium solution, we use the constant temperature com- 
mand to fix the temperature externally, allowing us to utilize 
Cloudy's sophisticated machinery to calculate cooling rates 
out of thermal equilibrium. In this manner, we create a grid 
of heating and cooling values as a function of temperature, 
gas density, chemical composition, and incident spectrum. 
The cooling rates presented in this work were created using 
version 07.02.01 of the Cloudy software. 

To automate the process of data production and orga- 
nization, we have created a code, called ROCO (Recursively 
Organized Cloudy Output.) ROCO uses a recursive algo- 
rithm to process user-specified loop parameters, making it 
possible to create data-grids of any dimension. Commands 
that are to be issued to Cloudy are given to the ROCO code 
in either one of two formats - loop commands with a set 
of parameters through which the code will iterate, and con- 
stant commands that are to be issued with the same value 
during each iteration over the loop commands. Since most 
uses of Cloudy involve the creation of large grids of models 
constructed by looping over a set of input parameters, the 
capabilities of ROCO give it the potential to be useful to 
a broader community of Cloudy users than just those who 
would use it to create the cooling tables discussed here. To 
this end, ROCO is structured in such a way that the post- 
Cloudy data analysis routines can be easily interchanged to 
suit the needs of different users. The code features an extra 
running mode that simply runs Cloudy over the specified 
parameter-space with no further processing of the data, as 
well as a template designed to help users create new run- 
ning modes suited to their specific needs. ROCO also has 
the ability to run multiple instances of Cloudy simultane- 
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ously, greatly reducing runtime. The parallel feature works iitot ~ uc/h- 
well on individual machines with multiple processors, as well 
as Beowulf clusters using the MPI (Message Passing Inter- 
face) framework. A copy of the ROCO code will be made 
available upon request to the authors. 

In Figure [T] we display the resulting cooling function for 
gas with uh = 1 cm"'' at metallicities, from Z = (metal- 
free) to 10 Zq, over the temperature range, 50 ^ T ^ 10* 
K. For these cooling rates, we use the coronal equilibrium 
command in Cloudy to simulate an environment free of ra- 
diation, where all ionization is coUisional. We also neglect 
the cooling from H2, so as to better illustrate the cooling 
contribution from metals at temperatures less than 10* K. 
We accomplish this by issuing the Cloudy command, no H2 H + e" - 

molecule. H" + H 



2.2 Implementation in Hydrodynamic Simulations 

We implement our metal cooling method in the Eule- 
rian adaptive mesh refinemen t hydrodynamic/N-b ody code, 
Enzo (|Brvan fc Normanlll997l : lO'Shea et al.ll2004 ). When a 
simulation is initialized, Enzo reads in the Cloudy/ROCO 
data-grid, storing the heating and cooling values as functions 
of temperature, H number density, and any other parame- 
ters, such as spectral intensity, depending on the nature of 
the simulation. The heating, F, and cooling, A, are stored 
with code units corresponding to [ergs s~^ cm^]. During the 
simulation, Enzo stores the mass density and internal energy 
for each grid cell in the box. At each hydrodynamic time- 
step, the radiative cooling solver cools the gas by lowering 
the internal energy via a simple Euler update, 

n + l n I • n , i-j. 



(1) 



where u"^^ denotes the internal energy of the grid cell with 
(x,y,z) coordinates, (i,j,k), at the (n-l-l)'th time-step, ii is the 
cooling rate in code units corresponding to [ergs s~^], and 
St is the adopted time-step. For every hydrodynamic time- 
step, the code subcycles through Equation [TJ selecting from 
three possible time-steps, until one full hydrodynamic time- 
step has been completed. The time-step, St, in Equation [T] 
adopts the minimum of the following three values: (1) half 
of the hydrodynamic time-step, (2) 10% of the cooling time, 
(u/u), or (3) the time remaining to have integrated over one 
full hydrodynamic time-step. The Euler update is the stan- 
dard method used for updating the internal energy in all of 
the established cooling routines in the Enzo code. Since the 
cooling rate is such a nonmonotonic function of temperature, 
this approach of using substeps with an explicit integration 
method has been fo und to yield the best com bination of 
accuracy and speed l|Anninos fc Normanlll993 ). The inter- 
nal energy and mass density for each grid cell are converted 
to temperature and number density and the heating and 
cooling values are calculated by linearly interpolating over 
the Cloudy/ROCO data-grid. The change in internal energy 
from the Cloudy/ROCO cooling rates is expressed as 



itc/R = (r - A)nH, 



(2) 



where uh is the H number density. 

We implement two distinct versions of the method de- 
scribed above. In the first and simplest version, the cooling is 
calculated solely from the Cloudy/ROCO data, as in Equa- 
tion (2] The total change in internal energy is 



(3) 

When converting the internal energy to temperature, it is 
necessary to know the value of the mean molecular weight, 
fi. In this implementation, we assume to be a constant 
with the value 1.22. For high temperatures, T > 10"* K, this 
method is sufficient for providing accurate gas cooling within 
hydrodynamic simulations. This implementation is not suit- 
able, however, when T < 10'' K and the formation of H2 be- 
comes important. Disregarding formation on grain surfaces 
and three-body formation, H2 primarily forms through the 
following channels: 
the H~ channel, 



H- + 7, 
> H2 e" 



and the HJ channel. 



H H"* 



+7, 



H+ -f H H2 + 



(4) 



(5) 



When a significant electron fraction exists, these reac- 
tions proceed to form H2 very quickly, with the H~ chan- 
nel typically dominating, except in the very high redshift 
universe ( z > 100), wher e H~ is readily de stroyed by 
the CMB llAbel et al.lll997l: iBromm et al.|[2002l ). Recently, 
iHirata fc Padmanabhan (|2006l ) have suggested that forma- 



tion of via 

He H+ ^ HeH+ + 7, 
HeH+ -I- H -» H+ -f He, 



(6) 



is responsible for more H2 than the HJ channel in Equa- 
tion [S] If, however, the ionization fraction is low, H2 forms 
very slowly, with equilibrium timescales that can exceed the 
current age of the universe. The consequence is that H2 for- 
mation is so sensitive to the thermal history of the gas that 
H2 fractions cannot be known without explicitly following 
the non-equilibrium chemistry during the simulation. We 
find this to be the case when using Cloudy to compute the 
cooling rate from H2. In searching for ionization equilibrium. 
Cloudy integrates over timescales that are unphysically long, 
leading to an overcalculation of the H2 fraction, producing 
cooling rates that are too high. 

We solve this problem in our second implementation by 
first removing the H2 molecule from Cloudy's chemical net- 
work with the no H2 molecule command. We, then, run two 
Cloudy/ROCO data-grids: one with the full set of elements 
and the other with H and He only. We obtain a metals-only 
data-grid by subtracting the metal-free data-grid from the 
compl ete data-grid. Using the establish ed H/He network in 
Enzo l|Anninos et all 1 19971 : 1 Abel et al.lil997) . we follow the 
non-quilibrium fractions of H, H+, H~, H2, Hj, He, He""", 
He"^^, a nd e~, and dir ectly calculate t he associated atomic 
(Black lS)8ll: ICen|[l99^ and molecular (|GaUi fc Pallalll998D 
cooling rate s. We provide the o ption to use the H2 cool- 
ing rates of iLepp fc Shulll (|l983t ). which are obsolete, but 
allow a means of comparison to older simulations. We also 
include the cooling, or heating, from electrons scattering off 
the CMB as 

Acomp = 5.4 X 10-^''(1 + zfn^{T - Tcmb), (7) 

where Tcmb = 2.7 (1 -f z) K l|Bromm et al.ll2002l ). We pre- 
vent the metals from cooling the gas below the CMB temper- 
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at ure by subtrac t ing th e metal cooling rate at T = Tcmb, as 
in lBromm et al.l l|200ll ). Including the CMB explicitly in cos- 
mological simulations would require adding an extra dimen- 
sion in the Cloudy/ROCO data-grid to account for the evo- 
lution of the CMB with redshift. While this is certainly pos- 
sible, interpolating over an extra dimension to calculate the 
cooling during the simulation would be significantly slower 
than the approximation described above. To test the va- 
lidity of this approximation of the CMB temperature-floor, 
we create a Cloudy/ROCO metals-only data-grid, explicitly 
including the CMB. At low densities (n ~ 1 cm~^), the val- 
ues of (A - r) from the data with the CMB included differ 
from the values of (A - A{Tcmb)) from the data without the 
CMB by a factor of roughly 2 near Tcmb- At higher tem- 
peratures and densities, the two values are nearly identical. 
The total rate of energy loss applied to the simulation gas 
in the second implementation is 



Utot = UH,He + UComp + / R, (8) 



where UH,He is the total atomic and molecular cooling from 
the H/He network, and uc/r is the metals-only cooling rate 
taken from the Cloudy/ROCO data in the manner described 
above. The value of fi is calculated directly from the H/He 
species fractions, neglecting any addition from the metals. 
In low-metallicity gases, this approach is reasonable, as the 
increase in jj. from the metals only reaches ~10~* for Z = 
10~^ Zq. Since the Cloudy/ROCO data-grids also store the 
values of for each point, this can be added to the value 
calculated without the heavy elements when the metallicity 
is very high. 

In Figure [S] we display low-temperature cooling func- 
tions for gases with varying density and metallicity, con- 
structed with the third implementation of the metal cooling 
method. To produce the data for Figure [21 we set up an un- 
physical, two-dimensional grid in Enzo that varies smoothly 
over density and temperature. We iterate the reaction net- 
work for a time equivalent to that between z = 99 and 
20 (~160 Myr), with hydrodynamics disabled, then com- 
pute the cooling with the third implementatio n of our metal 
cooling method, using the H2 cooling rates of ICalli fc Fallal 
l|l998h . Since the first stars are predicted to f orm at the cen- 
ters of ^10^ Mq dark matt er halos at z ~ 20 (|Tegmark et al.l 
I1997I : lYoshida et al.l l2003t ) , integrating the rate equations 
over this time interval places each of the species in the rel- 
ative abundances in which they would be found during the 
epoch of first-star formation. Hence, Figure [2] provides a di- 
rect comparison of the cooling rate of the gas that formed 
the first and successive generations of stars. The H2 frac- 
tions used to create the cooling rates for Figure [2] result 
from integrating the H/He chemical network for the period 
of time between z = 99 and 20. As such. Figure [2] should 
not be used as a general cooling function, except in the con- 
text mentioned above. However, Figures [3H8] do not include 
the cooling from H2 and may, therefore, be used as general 
purpose cooling functions when the cooling from H2 is not 
needed. 



3 METALS IN LOW-TEMPERATURE GASES 



3.1 Dominant Coolants 

Much attention has been given recently to the role of the 
first heavy elements in transitioning from the singular, 
high-mass mode of star formation of the first stars to the 
mode producing stars with a Salpeter initial mass fun ction 
(IMF)! Analytical stu dies by iBromm fc Loebl (|2003bl ) and 
ISantoro fc Shulll (|2006l ) have focused on the contributions of 
individual elem ents toward triggering f ragmentation in star- 
forming clouds. iBromm fc Loebl (l2003b) suggest C and O to 
be the dominant coolants in low-metallicity gas, in the pres- 
ence of an H2 dissociating UV b ackground created by the 
first stars (|Bromm fc Loebll2003al ). By calculating the cool- 
ing rate necessary to equate the cooling time to the free-fall 
time at n = 10* cm~^ and T = 200 K, the point where H2 
coolin g becomes ineffici e nt (lAbel et al.l I2OO2I : iBromm et al.l 
l2002h . iBromm fc Loebl (|2003bh predict individual critical 
abundances of C and O to be [C/H]crit — -3.5 and [0/H]crit 
~ -3.05, where [ A/HI = logio(NA/NH) - logio(NA/Nif)Q. 
ISantoro fc Shulll l|2006l ) consider the cooling from Fe and Si, 
in addition to C and O, and take into account the density 
dependence of metal cooling. In doing so, they find that 
the critical abundance of each element varies with density, 
reaching a minimum at a critical density that is different in 
each case. They also note that different elements dominate 
different density and temperature regimes. 

In Figs. |3Hni we plot the individual cooling contribu- 
tions for number density, uh = 10^ cm~'^ and metallicities, 
Z = 10~^ Zq and 1 Zq. In each case, we create a set of 
cooling data with the full complement of elemental species, 
from H through Zn, neglecting H2. We plot only the coolants 
whose contributions reach, at least, 10"'^ of the total cooling 
within the temperature range, 10 K ^ T ^ 5000 K. Since the 
cooling data were made assuming no incident ionizing radi- 
ation, as described in i|2.1l we observe the d ominant C tran- 
sitions to be from CI, instead of CI I, as in lBromm fc Loebl 
(|2003bl ) and ISantoro fc Shulll (|2006t ). For Z = 10"^ Zq, the 
most important coolants are the fine structure transtions at 
369.7 fim and 609.2 fim from CI and at 63.2 /xm from OI. At 
higher temperatures (T ^ 100 K), the Sill transition at 34.8 
^m becomes important as well. At higher metallicities, CO 
replaces atomic C (Fig.O and emission from Sil at 129.6 
becomes dominant (Fig.[5]|. Fe cooling is relatively unimpor- 
tant up to uh = 10^ cm^^, but is completely dominant for 
T ^ 200 K by uh = lO" cm"^, with OI strongest slightly 
below that temperature, and CO the most important be- 
low 50 K (Fig. E| ■ In addition to the elements covered by 
ISantoro fc Shuliri|2006l ). we note that cooling from neutral 
S reaches roughly the 10% level at uh = 10^ cm"^. The [SI] 
25.2 fj,m transition peaks at 40% of the total cooling at tin 
= 10^ cm"^ and T ~ 1000 K. The only other coolant that 
contributes at the level of at least 10"'^ of the total cooling is 
the 60.6 fj,m transition of [PII] at uh ~ W"' cm""'. The num- 
ber of coolants that reach 10"'' of the total grows quickly 
with density. We observe 23 distinct coolants contributing 
at that level at uh = 10^ cm"^, and 32 by 10^ cm"^. If we 
lower the threshhold to 10~^, there are a total of 84 coolants 
at nn ~ 10^ cm""', illustrating the strength of Cloudy and 
our cooling method. 
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3.2 Dust Grains 

Recently, studies have suggested that dust coohng at high 
densities can t rigger fragmentation for metaUicities as low 
as 10"^ ("Omukai et all l2005l: [Schneider et a lj l2006l: 

iTsuribe fc Om ukai 2006.: .Clark et al.l l2007'). Schnei der et al l 
have claimed that between 15 and 30% of the mass 
of the progenitor of a pair-instability supernova is con- 
verted in dust. Ho wever, observations of the Crab nebula 
ijCre en ct al.' 2004) and the Cassiopeia A supernova rem- 
nant (Krausc ct al.. ,2004 ) have returned little or no signs of 
dust, suggesting that type II supernova may not actually 
be large dust producers. Given the controversy surrounding 
the existence of dust grains in the formation environments of 
second-generation stars, we do not include the cooling from 
dust in the analysis, but leave it for a separate work. To 
provide an example of the effect of dust on the cooling rate, 
we run a simple model including dust in Cloudy. We use a 
model designed to simulate the dust within the ISM, using 
the Cloudy command, grains I SM. The dust ph ysics used in 
Cloudy is described in detail bv lvan Hoof et al.i ( ,2004 ) ■ The 
ISM dust model in Cloudy consists of both graphite and sil- 
icates with sizes ranging from 5 x 10~^ fim to 0.25 /im and 
a power-law size-distribution with a power-law index of -3.5. 
For solar metallicity, the total grain abundances per H are 
■[^g-9.811 £qj. gj-apijjtg and lO"*'-^"'* for silicates. In Fig.[8l we 
plot the cooling rate from metals at uh ~ 10^ cm""^, with 
and without dust grains, for metaUicities, Z = 10~^ Z©, 10~* 
Zq, and 10~^ Zq. Since it is inappropriate to think of the 
dust and gas-phase metal abundances as independent, we 
directly scale the dust abundances with the gas-phase metal 
abundancs. For number densities lower than 10® cm~^, the 
additional cooling from dust is almost negligible. At higher 
densities, dust becomes more important. If dust is, in fact, 
produced in the supernovae of the first stars, it is likely to 
be as important as previous studies have claimed. 



3.3 Thermal Instability and Fragmentation 

In order to study the ability of a collapsing gas-cloud to 
fragment, we first identify the regions of density and tem- 
perat ure where t he classical fragmentation criterion, tcooi < 
tdyn l|Fieldlll965h . is met, with the dynamical time expressed 



tdyn 



16Gp' 



(9) 



We limit this analysis to solar abundance patterns. This rep- 
resents the first step of an incremental approach to studying 
the general criteria that lead to fragmentation in collapsing 
clouds. The use of solar abundance patterns will allow us 
to begin to quantify the chemical abundance required for 
fragmentation. In a following paper, we will study nonsolar 
abundance patterns motivated by predicted yields of pri- 
mordial supernovae. In this future work, we will vary the 
abundances of individual elements, which will require the 
exploration of a much larger parameter-space. We create 
a Cloudy/ROCO data-grid with the following parameters: 
50 K < T < 1000 K with 51og(T) ~ 0.012 (100 points), 1 
cm"^ s: n^r < lO" cm"^ with 5log{nH) = 0.1, and 10"'' 
Zq Z ^ 10"^ Z© with 51og(Z) = 1. We, then, follow 
the same procedure used to produce Fig. (2] as described 



in i]2.2l In this section, we omit the decrease in the cooling 
rate caused by Compton heating on the CMB. In Fig. |9] we 
plot the log of the ratio of the dynamical time to the cool- 
ing time for each of the metaUicities in the data-grid and 
for the metal- free case. A cloud is able to fragment when 
logio(ttjy„/tcooi) > 0. As expected, there is no density at 
which metal- free gas can fragment for T < 200 K. As the 
metallicity increases, the value of \ogio {tdyn /tcooi) slowly in- 
creases, first in the low-temperature regime, where H2 cool- 
ing is inefficient, so even a small amount of metals has an 
effect. For gas with uh = cm"^ at T = 200 K, the frag- 
mentation criterion is nearly met by Z = 10~® Zq and well 
satisfied one order of magnitude higher in metallicity. Once 
the metallicity reaches 10~^ Zq, the entire parameter-space 
is fragmentable. At high densities, however, fragmentation 
will be curtailed as the clou d becomes optically thick to 
its own radiation iILow fc Lvnd cn-BcU 1976; Rccs 1976). As 
was also reported bv lSantoro fc ShuUl (|2006i ). the efficiency 
of the metal cooling peaks at riH ~ 10^ cm~^, significantly 
lowering the critical metallicity required for fragmentation. 

The addition of metals to a gas also has the potential 
to trigger therma l i nstab i lities during cloud-collapse. As in 
lAbel et all (|2002l ) llieii [l96i) , we define a parcel of gas 
losing energy at a rate, L, to be thermally unstable if 



+ L(p,T) > 0, 



(10) 



where L is expressed in terms of the cooling rate. A, as 

L(p,T) =pA(p,T). (11) 

The cooling rate. A, can be locally approximated by a power- 
law in both temperature and density as 

The partial derivatives of Eon. llOl become 



(/3 + i)A(p,r) 



and 
(-) 



T 



The thermal instability criterion simplifies to 
Q - /3 < 2. 



(13) 
(14) 
(15) 



In Fig. [To] we plot the value of the instability parameter, 
(a - /3), for the same cooling data used for Fig. |9] For 
metal- free gas (Fig. 1101 top-left), the instability parameter is 
greater than 3 over nearly the entire par ameter spa c e, and 
alw ays greater than 4 at high densities. lAbel et al.l (|2002l ) 
and IVoshida et all (I2006D both arrive at the same conclu- 
sion using this analysis for metal- free gas, with the added 
assumption that the cooling function is independent of den- 
sity. As the metallicity increases, a region of thermal in- 
stability forms at low density and temperature. When the 
metallicity reaches 10~* Zq (Fig. 1101 middle-right), a sec- 
ond unstable region exists for 10^ cm~'^ < nn < 10^ cm~"^, 
at a temperature of a few hundred K. The second unstable 
region coincides with the increase in cooling efficiency to its 
maximum value, illustrated in Fig. [51 

Fragmentation is more likely to occur when the gas is 
both thermally unstable, and can cool faster than the dy- 
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namical time. We indicate the regions where both the frag- 
mentation and thermal instability criteria are met in white 
in Fig. 111! No doubly unstable realm exists for metallicities, 
Z ^ 10"^ Zq. 



3.4 Effects of the Cosmic Microwave Background 

The CMB creates a temperature floor, below which gas can- 
not cool radiatively. We study the influence of the CMB on 
the evolution of star-forming gas by applying a CMB floor at 
z = 20 in the manner described in H2.2l to the cooling data- 
grid used in i|3.3l The CMB affects the cooling properties of 
the gas in two ways. The first is by increasing the cooling 
time at temperatures near the CMB temperature to greater 
than the dynamical time so that the fragmentation criterion 
is no longer satisfied. The second is by increasing the value 
of a, from Equation 1151 at low temperatures, making the 
gas thermally stable. In Fig. 1121 we illustrate the influence 
of the CMB on the doubly unstable regions, shown previ- 
ously. The unstable region that existed in the low-density, 
low-temperature regime is completely eliminated. There re- 
mains, however, a small area of instability for metallicities 



4 DISCUSSION 

We have introduced a new method for including the ra- 
diative cooling from metals in large, three-dimensional hy- 
drodynamic simulations. In addition to its implementation 
in the AMR code, Enzo, t his method has also been used 
bv iBogdanovic et all l|200d ) in numerical simulations with 
the smoothed particl e hydrodynami cs (SPH) code. Gadget 
(jSprineel et al.ll200ll : ISpringeJ lioosh . Our technique takes 
advantage of the extremely complex chemical reaction net- 
work of the preexisting radiative transfer code. Cloudy, 
which includes a full elemental coverage from H to Zn, along 
with a variety of molecular species and dust. With the sin- 
gular assumption of ionization equilibrium for the heavy el- 
ements, we are able to precalculate cooling rates for gases 
with any chemical abundance in all manners of radiation en- 
vironments over a temperature range of 10 to 10* K. With 
cooling rates computed in advance, we eliminate the bar- 
rier that has classically prevented large chemical models 
from being incorporated into three-dimensional numerical 
simulations. Because our cooling scheme is valid over such 
a large range of density and temperature, and features so 
many coolants, it can be applied to a huge variety of as- 
trophysical problems, such as the evolution of the ISM and 
IGM, normal star formation, planetary nebulae, accretion 
disks, and protoplanetary disks. 

One advantage of the large chemical network of Cloudy 
is that we are able to determine the dominant coolants from 
a complete sample of atomic species up to an atomic number 
of 30. Fine-structure transitions of C and O are the greatest 
contributors to the cooling up to number densities of about 
10^ cm~^, where Fe cooling becomes significant and C is 
marginalized. The importance of these three elements, along 
with Si, in triggering the formation of t he first low-mass 
stars has been studied in great detail by ISantoro fc Shulll 
(|2006l ). The cooling models we have constructed using so- 
lar abundance patterns reveal S cooling to be important for 



cm ^. S is produced on ly slightly 



10^ cm-^' < riH < W . 

less than Si in type II llWooslev fc Weaver 19951 ) and pair- 
instability supernovae (|Heger fc Wooslevll2002i r and should 
be taken into account when considering the metals respon- 
sible for the transition from Population III to Population II 
star formation. The ability to specify individual abundances 
in our cooling method makes it straightforward to simulate 
the evolution of gas with non-solar abundance patterns. 

In addition to elevating the cooling rate of a gas to 
satisfy the classical fragmentation criterion, metals also in- 
crease the potential for fragmentation by creating thermal 
instabilities. We have identified regions in temperature and 
density in which both the classical fragmentation and ther- 
mal instability criteria are met to be the physical conditions 
most likely to see fragmentation occur. We observe these 
doubly unstable regions to exist for metallicities as low as 
10^* Zq. If fragmentation cannot occur outside these re- 
gions, then the fate of a star-forming gas-cloud will be deter- 
mined by the path taken through density-temperature space 
as it collapses. If we consider the doubly unstable regions 
in Fig. 111! appropriate for star-formation in current epoch, 
there will almost certainly be a period of double instability 
when Z ^ 10"'' Zq. Int erestingly, the densit y-temperature 
tracks shown in Fig. 1 of lOmukai et al.l |200^ indicate that 
gas with Z — 10"'* Zq will pass rig ht through the stable 
region that separates two instabilities. lOmukai et al.l (|2005l ) 
also find that star formation at that metallicity only pro- 
duces high-mass fragments. At high redshift, the CMB sig- 
nificantly reduces the size of the doubly unstable regions. 
Thermal instabilities, though, are extremely sensitive to the 
slope of the cooling rate as a function of density and tem- 
perature. Every element has distinct cooling properties, and 
will, therefore, produce different thermal instabilities. As 
such, the key to uncovering the nature of the first Population 
II stars will be in the determination of the mass function of 
their precursors. 



5 FUTURE DEVELOPMENT 

In papers to follow, we will extend our study to gases with 
non-solar abundance patterns. We will explore thermal and 
double instabilities created by individual elements, as well as 
abundance patterns produced in Population III supernovae. 
Future studies will also investigate the effects of background 
radiation on the evolution of star-forming gas. The final 
word, however, will only come from three-dimension al, nu- 
merica l simulations. The simulations by Smith & Sigu rdssonI 
(j2007l ). employing the methods described here with solar 
abundance patterns, have confirmed that fragmentation oc- 
curs for metallicities, Z ^ 10"^ Zq. We intend to pair all 
future predictions made from analysis of thermal instabili- 
ties with full numerical simulations. 

Although dust physics has been implemented in Cloudy, 
we currently treat only metals in the gas-phase in our anal- 
ysis. In the future, we will study the effects of dust cooling 
in more detail. We will also couple the dust chemistry to the 
H/He chemical network in Enzo, so as to properly model the 
formation of II2 on grain surfaces. A great strength of this 
method is the use of the ever-expanding chemical network of 
Cloudy. As more physical processes are incorporated into the 
Cloudy software, the utility of this method will increase as 
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well. One major constraint of this work, however, is that its 
validity is confined to the optically thin limit. The approx- 
imations made here break down at opacities of order unity. 
For higher opacities, our method provides a core module 
for flux-limited diffusion schemes. Complex geometries in- 
troduce problems of self-heating and shadowing, which will 
require full, three-dimensional radiative transfer. 
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Figure 1. Cooling functions, excluding cooling from H2, for gases 
with iiH = 1 cm^'^ and metallicities, Z = (solid), 10""^ Zq (dot), 
V)-'^ Zq (dash), 10"^ Zq (dot-dash), 1 Zq (dot-dot-dash), and 
10 Zq (dot-dash-dash). 
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Figure 2. Cooling functions, including H2 cooling, for gases with 
n^f = 1 cm^"^ (top-left), iiH = 10'^ cm^'^ (top-right), iih = 10® 
cm""^ (bottom-left), and nj/ = lO" cm^'^ (bottom-right). Metal- 
licities are Z = (solid), IQ-'^ Z© (dot), lO'^ Zq (dash), IQ-* Z© 
(long dash), lO'^ Z© (dot-dash), and lO'^ Z© (dot-long dash). 




Figure 3. Cooling contributions from C and O species that reach 
at least 10^"^ of the total cooling for gas with with iij/ = 10'"'' 
cm~^ and Z = 10"'^ Zq. The total cooling (solid, black) includes 
all species contained within the Cloudy chemical network. Com- 
ponents shown are [Ol] 145.5 fim (dot-dot-dash), [01] 63.2 
(dash-dash-dot), [CI] 369.7 ^tm (dash), [CI] 609.2 (long dash), 
CI 985 nm (dash-dot), and [CII] 157.6 f^m (long dash-dot). The 
solid, grey line indicates the sum of all the components plotted. 




Figure 4. All other coolants not plotted in Fig. |3] that reach 
at least 10~^ of the total cooling for gas with njj = lO'^ cm~^ 
and Z = 10"'^ Zq. The total cooling (solid, black) includes all 
species contained within the Cloudy chemical network. Compo- 
nents shown are [Fell] (dot), [SI] 25.2 fim (dash), [SI] 56.3 fiia 
(long dash), [Sil] 129.6 ^lm (dash-dot), [Sil] 68.4 ^m (long dash- 
dot), and [Sill] 34.8 fim (dash-dot-dot). The solid, grey line indi- 
cates the sum of all the components plotted. 




Figure 5. Cooling contributions from C and O species that reach 
at least 10^"^ of the total cooling for gas with with iij/ = 10'^ 
cm~^ and Z = 1 Zq. The total cooling (solid, black) includes 
all species contained within the Cloudy chemical network. Com- 
ponents shown are CO (dot), [OI] 145.5 ^m (dot-dot-dash), [01] 
63.2 (dash-dash-dot), [CI] 369.7 ^lra (dash), [CI] 609.2 
(long dash), and CI 985 nm (dash-dot). The higher dotted line 
represents CO emission from C^'^O^'', while the lower dotted line 
shows emission from C^'^O^*'. The solid, grey line indicates the 
sum of all the components plotted. 




Figure 6. All other coolants not plotted in Fig. |5] that reach at 
least 10~^ of the total cooling for gas with n^f = 10'^ cm~^ and Z 
= 1 Zq. The total cooling (solid, black) includes all species con- 
tained within the Cloudy chemical network. Components shown 
are [Fell] (dot), [Fel] 24.0 (dash), [Fel] 34.7 /^m (long dash), 
[SI] 25.2 /im (dash-dot), [SI] 56.3 ura (long dash-dot), [Sil] 129.6 
lira (dash-dot-dot), and [Sil] 68.4 /im (dash-dash-dot). The solid, 
grey line indicates the sum of all the components plotted. 
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Figure 7. Subset of the most dominant coolants at n// = lO'* 
cm~^ and Z = 1 Z©. The total cooling (solid) includes all 
species contained within the Cloudy chemical network. Compo- 
nents shown are CO (dot), [Fel] 24.0 fim (dash), [Fell] (long 
dash), [OI] 63.2 fj,m (dash-dot), and OI 630 nm. 
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Figure 8. Total cooling rate from metals for gas at nu = 10 
cm^'^, with metallicities Z = 10~^ Zq (solid), 10~^ Zq (dashed), 
and 10~s Zq (dash-dot). The black lines indicate the total cool- 
ing from gas-phase metals only. The grey lines show the cooling 
with gas-phase metals and dust grains, created with the ISM dust 
grain model in Cloudy, using the command, grains ISM. The ISM 
dust model in Cloudy consists of both graphite and silicates with 
sizes ranging from 5 X 10^'^ fim to 0.25 fim and a power-law size- 
distribution with a power-law index of -3.5. For solar metallieity, 
the total grain abundances per H are IQ^"-**!! for graphite and 
-|^Q— 9.748 silicates. In each case shown, the dust grain abun- 
dances have been scaled to the gas-phase abundances. 
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Figure 9. Contours of logio (t|jy„/tcooi) over number density and 
temperature for gases with metallicities, Z = (top-left), 10^® 
Zq (top-right), 10-5 Zq (middle-left), lO"-* Z© (middle-right), 
10"^ Zq (bottom- left), and 10"^ Zq (bottom-right). H2 cooUng 
is included. 
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Figure 10. Contours of the instability parameter, (q - /3) over 
number density and temperature. The medium is unstable for 
values less than 2. The metallicities are Z = (top-left), 10~® Zq 
(top-right), 10-5 Zq (middle-left), 10"" Zq (middle-right), 10-^ 
Zq (bottom-left), and 10"^ Z© (bottom-right). At Z = lO""* Z©, 
two separate thermally unstable regions exist. These two regions 
merge by 10~^ Zq. 
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Figure 11. The white patches indicate regimes of density and 
temperature where logiQ{tj^y„/t^^^i) > and {a - (3) < 2 for 
metallicities, Z = lO"* (top), 10~^ Zq (middle), and lO'^ Zq 
(bottom). As in Fig. 1101 there are two individual doubly unstable 
regions at Z = 10~^ Zq that have merged by 10"'^ Zq. 
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Figure 12. Doubly unstable regions for the same metallicities as 
in Fig. 111! but with a CMB temperature floor at z = 20 included. 



